#!/bin/env Rscript

library(HMMcopy)

argss <- commandArgs(trailingOnly = TRUE)

print(argss)




samplew <- argss[2]
controlw <- argss[3]
gcw <- argss[4]
mapw <- argss[5]

# Correct readcount and estimate copy number
sample_uncorrected_reads <- wigsToRangedData(samplew, gcw, mapw)
control_uncorrected_reads <- wigsToRangedData(controlw, gcw, mapw)
sample_corrected_copy <- correctReadcount(sample_uncorrected_reads)
control_corrected_copy <- correctReadcount(control_uncorrected_reads)

# Segmentation

# Normalizing sample by control
sample_corrected_copy$copy <- sample_corrected_copy$copy - control_corrected_copy$copy
# Segmenting
segmented_copy <- HMMsegment(sample_corrected_copy)
# retrieve converged parameters via EM
# param <- HMMsegment(sample_corrected_copy, getparam = TRUE)
# param$mu <- log(c(1, 1.4, 2, 2.7, 3, 4.5) / 2, 2)
# param$m <- param$mu
# perform segmentation via Viterbi
# segmented_copy <- HMMsegment(sample_corrected_copy, param)


# Export

# Export to SEG format for CNAseq segmentation
# rangedDataToSeg(sample_corrected_copy, file = "sample_corrected_copy.seg")

# Visualization

chr <- c(1:22, "X", "Y")

# pdf("1.hmmcopy.pdf")

# plotBias(sample_corrected_copy)
# plotCorrection(sample_corrected_copy)
# plotSegments(sample_corrected_copy, segmented_copy, lwd = 4)
source(paste(Sys.getenv('tools_path'), 'script/cnv_p_hmmcopy.r', sep='/'))
cnvphmmcopy(sample_corrected_copy, segmented_copy, chr=chr)

# dev.off()

q(save = "yes", status = 0, runLast = FALSE)
